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Abstract 

Recent theoretical investigations on quantum coherent feedback networks have 
found that with the same pump power, the Einstein-Podolski-Rosen (EPR)-like en¬ 
tanglement generated via a dual nondegenerate optical parametric amplifier (NOPA) 
system placed in a certain coherent feedback loop is stronger than the EPR-like en¬ 
tangled pairs produced by a single NOPA. In this paper, we present a linear quantum 
system consisting of two NOPAs and a static linear passive network of optical de¬ 
vices. The network has six inputs and six outputs, among which four outputs and 
four inputs are connected in a coherent feedback loop with the two NOPAs. This 
passive network is represented by a 6 x 6 complex unitary matrix. A modified steep¬ 
est descent method is used to find a passive complex unitary matrix at which the 
entanglement of this dual-NOPA network is locally maximized. Here we choose the 
matrix corresponding to a dual-NOPA coherent feedback network from our previous 
work as a starting point for the modified steepest descent algorithm. By decomposing 
the unitary matrix obtained by the algorithm as the product of so-called two-level 
unitary matrices, we find an optimized configuration in which the complex matrix is 
realized by a static optical network made of beam splitters. 


1 Introduction 

In recent years, research related to the Einstein-Podolski-Rosen (EPR) entanglement in 
continuous-variable quantum information processing has become increasingly vital since it 
can be shared by two distant communicating parties and is used as the crucial resource for 
important applications such as quantum teleportation and superdense coding [I]. Com¬ 
pared to discrete variable entangled states, continuous variable entanglement such as EPR 
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entanglement [T] is generated efficiently by a pair of squeezed light beams and utilized 
with expeditiousness in measurement of quantum states which is a critical step in quan¬ 
tum communication protocols M. 

A device that is used to produce EPR-like entangled states is a nondegenerate opti¬ 
cal parametric amplifier (NOPA), which contains a cavity with a nonlinear crystal 
inside. Via a strong undepleted coherent beam pumped to the crystal, two ingoing sig¬ 
nals in vacuum state interact with two modes of the cavity separately, and generate two 
output beams which are squeezed in quadrature-phase amplitudes and considered as EPR 
entanglement [lj. EPR entanglement between the two outgoing fields is measured by the 
two-mode squeezing spectra of the fields. A strong EPR entanglement is denoted by a high 
degree of two-mode squeezing. As a quantum system in reality is sensitive to its external 
environment, it undergoes unwanted interaction with an external electromagnetic held as 
a “heat bath” [5]. Therefore the strength of EPR entanglement generated by such an open 
system can be degraded due to transmission losses and decoherence, which leads to limited 
communication distance. Thus methods to enhance EPR entanglement are of interest to 
improve the quality of quantum communication. 



Figure 1: The dual-NOPA coherent feedback network. 


fin, 1 fin,2 



fout,l font, 2 

Figure 2: System structure. 


Coherent feedback is a feedback loop which directly connects quantum optical compo¬ 
nents without employing any measurement apparatus in the loop [B]. Our previous work 
[7] presents a dual-NOPA coherent feedback scheme comprised by two NOPAs as shown in 
Fig.0 Each NOPA Gi (i = 1,2) in the figure is simply denoted by a block with four inputs 
and two outputs. Gi contains two modes a* and bj. Ingoing signal and amplification 
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loss £,ioss,a,i interact with mode a*; similarly and amplification loss £,i OS s,b,i interact with 
mode bi . The output £out,a, 1 of Gi is connected to the input ^ ni0j2 to G 2 , and the outgoing 
signal C, outt b ,2 from G 2 is the input to G\. EPR entanglement is generated between outputs 
iout,a ,2 and £out,b,i- A more detailed description of the NOPA will be given in Section 2.3 


EPR entanglement generated from a dual-NOPA coherent feedback network where two 
NOPAs are placed at two endpoints (Alice and Bob) separately is compared to that of a 
single NOPA located in the middle (at Charlie’s), at a location between Alice’s and Bob’s, 
in [?]. When amplification losses are neglected, under the same values of configuration pa¬ 
rameters such as transmission losses, decay rates and pump power, the coherent feedback 
network improves EPR entanglement between two outgoing fields in terms of increasing 
the level of two-mode squeezing that can be achieved over the single NOPA, see [7j. Also, 
with the same setting of decay rates and when losses are neglected, the coherent feedback 
network requires less pump power to generate the same degree of two-mode squeezing com¬ 
pared to the single NOPA. Thus the coherent feedback scheme has improvement in EPR 
entanglement generation. 


Based on the above facts, we consider a system consisting of two NOPAs connected in 
a coherent feedback loop with a static passive linear network which is denoted by a 6 x 6 
complex unitary matrix S, as shown in Fig. [ 2 ] The passive linear network can be assembled 
by several static linear optical devices, such as beamsplitters and phase shifters, of which 
the transformation functions are 2x2 unitary matrices [H]. The EPR entanglement 
between the continuous-mode fields £ 0 ut.i and £ 0 «t ,2 is of interest. The coherent feedback 
network in Fig. [I] corresponds to a special case where S takes on a particular value that 
will be given in Eq. 0 in Section [3] EPR entanglement is quantified by the amount of 
two-mode squeezing between the two fields at the frequency u = 0 rad/s. The two-mode 
squeezing will be given by a certain nonnegative-valued function V (0; S) of the matrix 
S (to be defined in Section II-B), and strong EPR entanglement corresponds to a small 
value of this function. Thus, the aim of the paper is to optimize the two-mode squeezing 
by finding a local minimum (denoted by 5) m ) of the real-valued function V(0; S) subject 
to the constraint that S is unitary, which can be solved by a modified steepest descent 
algorithm on a Stiefel manifold as proposed in [15], with S c fb as an initial point. Via the 
decomposition of 5) m into the product of 15 two-level unitary matrices [IS], we can then 
find the configurations of optical devices that realizes the passive linear network S) m . 

The structure of the rest of this paper is as follows. A brief review about linear quantum 
systems, EPR entanglement between two continuous-mode fields and dynamics of a NOPA 
is given in Section [2} Section [3] describes the system of interest and Section [4] explains 
the optimization process. In Section [BJ by a decomposition of the unitary matrix S) m , a 
detailed physical configuration of the whole network is presented. Finally Section [6] gives 
a short conclusion of this paper. 
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2 Preliminaries 


The notations used in this paper are as follows: i = \/—1 and Re denotes the real part of 
a complex number. The conjugate of a matrix is denoted by ■#, - T denotes the transpose 
of a matrix of numbers or operators and •* denotes (i) the complex conjugate of a number, 
(ii) the conjugate transpose of a matrix, as well as (iii) the adjoint of an operator. O mxn 
is an m by n zero matrix and I n is an n by n identity matrix. Trace operator is written 
as Tr[-] and tensor product is <g). Also, eig(-) denotes eigenvalues of a matrix and max(-) 
denotes the maximum value. 

2.1 Linear quantum systems 

An open linear quantum system without a scattering process contains n-bosonic modes 
ctj(t) (j — 1,..., n) satisfying [oj(t), aj(t)*] = Sij. The system interacts its environment via 
a time-varying interaction Hamitonian 


m 


Hint(t) = iy - Lfcw), 


3 =1 


(1) 


where Lj is the j-th system coupling operator and Zj{t) (j — 1,..., m) is the field operator 
describing the j'-th environment field [5]. When the environment is under the condition of 
the Markov limit, the field operators satisfy [£j(t), ^-(s)*] = S(t — s), where S(t) denotes the 
Dirac delta function. When Lj is linear and H is quadratic in cq and a*, the Heisenberg 
evolutions of mode a 3 and output filed operator ^ out: j are defined by cij (t) = U(t)*a,jU(t) 
and iout.j{t) = U(t)*£i n> j(t)U(t ) with unitary U(t) = ex^ (—i H int (s)ds) and satisfy [8J, 


where 


with quadratures mi 


z 

z 

£,out 




z(t) = Az(t) + B£(t), 
Zout,j(t ) = Cz(t) + D£(t). 


(af,af,...,a«,<) T , 

(f q F p f q f p ) T 

\^out, 1> Sont.l; ■ ■ ■ i s out,li ’ out,U ’ 


aj + a*, a p = (aj - a*)/i, 

Zj+Zj, Z’j (s, -s?)/h 


( 2 ) 

(3) 


(4) 


(5) 


2.2 EPR entanglement between two continuous-mode fields 

It is important to note that the output fields £^1 and Zout ,2 are two continuous-mode 
Gaussian fields rather than two single mode Gaussian fields. That is, each of Zout, 1 and 
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£out ,2 contain a continuum of modes rather than just a single mode. Therefore, the en¬ 
tanglement of the fields cannot be measured using entanglement measures for bipartite 
Gaussian systems, such as the well-known logarithmic negativity measure ms- Instead, 
when the incoming fields are in the vacuum state, the EPR entanglement of f out ,i and Gut ,2 
is assessed in the frequency domain using two functions V+{iuj) and V-{ioj ) cam nu that 
will be detailed below . 

F(ioj), the Fourier transform of f(t ) is defined as F (iu) = JZ f (t) e~ VjJt dt. Sim¬ 
ilarly, we get the Fourier transforms of £ 0 «ui(f), Cut 2(^)5 z(t) and £(t), as S outj i (xui), 
Eout ,2 Z(iuj) and S (iu), respectively. Using 0 and the definition of the Fourier 
transform, we get 

S L,lM + = 1-0 o CLt^e^dt + I-oo Cut, 2 ( 1 )^ dt = C ± Z M + D l E M , 

- S L, 2 M = IZ Cut, i(t) e ~ tujt dt - = c 2 z (tu) + u 2 “ m ,(6) 

where Ci = [1 0 1 0]C, C 2 = [0 1 0 —1 ]C, D\ = [1 0 1 0 ]D and D 2 = [0 1 0 —1 ]D. 

The two-mode squeezing spectra V + (iu) and V_{ioj) are defined as 

ftiH + + SL, 2 K))) = - "'), 

fcH - - ^, 2 K))> = U-M<^ - u/), (7) 

where (•) denotes quantum expectation. V + [iuj ) and U_(za;) are real valued and can be 
easily calculated as described in Haim 

U+M =Tr [H^iuj)*H^iu )], 

U_M =Tr[H 2 (ioj)*H 2 (iuj)\, 

where H\ and H 2 are transfer functions 

Hj(uj) = Cj {iul — A)~ l B + Dj, (j = 1,2). (10) 

Denote V(itu) = V + {ioj ) + V_{ioj). A sufficient condition for the fields Cut,\ and Cut ,2 to be 
EPR-entangled at the frequency u> rad/s is [TT], 

V(iu) < 4. (11) 

Ideally, we would like V(iuj) = V±{mf) = 0 for all u>, which denotes infinite-bandwidth two¬ 
mode squeezing, representing an ideal Einstein-Podolski-Rosen state. However, in reality 
the ideal EPR correlation can not be achieved, so in practice the goal is to make V(iuf) as 
small as possible over a wide frequency range HU- 

Define Cut,i = ^Cout, 1 , Cut ,2 = F^ 2 Cut ,2 with -0i,-0 2 G (—7T, 7 t] and denote the corre¬ 
sponding two-mode squeezing spectra as V± 1 , ^ 2 ), we have the following definition 

of EPR entanglement. 

Definition 1 Fields Cut, 1 and Cut ,2 are EPR entangled at the frequency u rad/s if 3 ifi, if 2 6 
(— 7 r, 7r] such that 

(ilu, ip \, ip2 ) T {tat, ip \, iff) <4. (12) 


( 8 ) 

(9) 
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Unless otherwise specified, throughout the paper, EPR entanglement refers to the case with 
ifi — '02 = 0 . EPR entanglement is said to vanish at u ; if there are no values of if \ and if 2 
satisfying the above criterion. 

Following BE], we have a good approximation V + {iuj) ~ ld + (0) and V-(iu>) ~ Vh(0) 
at low frequencies. Thus in the rest of paper, we focus on V±(0) at frequency u = 0 as a 
measure to quantify EPR entanglement. 

2.3 The nondegenerate optical parametric amplifier (NOPA) and 
the dual-NOPA coherent-feedback network 

A NOPA ( Gi ) is a linear quantum system with four ingoing fields in the vacuum state and 
two outgoing fields, as shown in Fig. [ 3 } By shining the pump beam onto the x ^ nonlinear 


/. ... 




NOPA 


A . Uncii 

7 loss,a,i 

Gi 


r fiKb, 

^ loss,b,i 

a i b i 



Figure 3: Input/output block representation of a NOPA. 


crystal, the two oscillator modes a, and bi inside the cavity become coupled via the two¬ 
mode squeezing Hamiltonian H = |e(a*fe* — a A), where e is a real coupling coefficient 
related to the amplitude of the pump beam [4]. The modes satisfy the commutation 
relations [a*, a*] = 1, [bi,b*] = 1, [a*, 5.;] = 0, and [a;,&*] = 0 |5j. Mode is coupled to 
ingoing noise £,i n , a ,i and amplification loss fioss,a,i via the coupling operators L\ = 
and L 3 = y Ufa, respectively, for some constant damping rates 7 and At; similarly mode bj 
interacts with input signal and additional noise fi oss ,b,i by operators L 2 = y/fib and 
L 4 = \fkb. The dynamics of the NOPA ( Gi ) is 


di ( t ) 


bi (t) 


^ r, ^ a i (f) T (f) \ffif,in,a,i if) V~iff,loss,a,i (f) , 

^ bi ft) T “Uj (b) ft) \fnf,loss,b,i ft) , 


7 + K 


(13) 


with outputs 


f,out,a,i (f) —yfiyCH (t) T f,in,a,i ft) , 

£out,b,i ft ) =\ffbi ( t ) + £,in,b,i if) ■ (14) 


More details of the standard NOPA set up can be found in [4]. 
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3 The system model 


We consider the system shown in Fig. [2} The whole network consists of two NOPAs and a 
static passive linear sub-system denoted by S. The sub-system has six input signals £m,i, 
fin,2, fout.a. 1 , £out,b, i, font,a,2 and £ ou t,b, 2 , among which Cout,a,i, £out,b ,i are outputs of NOPAl, 
and font,a, 2 , font, 6,2 are outputs of NOPA2. Among the six outgoing beams of the passive 
sub-system are the two ingoing signals into NOPAl £ in>ai i and the two input beams 

of NOPA2 £i nt a ,2 and fm,b, 2 , and the EPR entanglement of interest is between f^p and 
font, 2 - The ingoing fields of the passive linear network are all in the vacuum state [13]. The 
transfer function of the static sub-system is a 6 x 6 matrix S, thus we have 


(15) 


£,out,l 


Cm,l 

£,out, 2 


Cm,2 

Cm,a, 1 
Cm,6,1 

= s 

Cowt,a,l 

Coixt,6,l 

Cm,a,2 


£,out,a,2 

Cm,6,2 


£,out,b, 2 


where S is a complex unitary matrix m, 

S*S = SS* = I 6 . 


(16) 


For the static passive linear matrix for the dual-NOPA coherent feedback network [7] 
shown in Fig. [lj the matrix S is 


S c fb - 


0 

0 

1 

0 

0 

0 


(17) 


Both NOPAs (Gi and G 2 ) in the network are identical with the same coupling constants 


e, 7 and k as discussed in Section [2~3j G 1 has two modes eq and b±, and G 2 contains modes 
a 2 and b 2 ■ The oscillation modes follow the commutation relations [a*, a*] = 5 i3 , [a*, b 3 ] = 0, 
[a*, b*\ = 0, [aj, a,j\ = 0 and [b u b 3 ] = 0 (i,j = 1 , 2 ). 

Define the following quadratures 

z =K, a p 3 , b\, &£, a q 2 , a p , b\ , b^f, 

f = r ci cp tP fP tP tg tP 1 t 

<,loss l^>loss,a,l 5 ^ loss,a,1 ’ ^loss, 6,1 ’ ^loss, 6,1 ’ ^ loss,a,2 ’ ^ loss,a,2 5 ^ loss,b,2 ’ ^loss,b,2) 5 

AO = rt<? tP tg tP ]T 

> LS>m,l ? Sin,1 ’ >m,2’ >m,2J ’ 

d(°) _[><? dP dP dP lT 

f=[f (i)T ,f£j T - (18) 

Define the real unitary matrix S' as the quadrature form of matrix S'. Based on the defining 
equations for the quadratures ((Hj), the relation between S and S is 

(19) 


S = ^ KSK* + ^K*S*K t , 
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where 


K = I 6 <g> 


1 

—i 


( 20 ) 


According to the dynamics of the two NOPAs given by (13) and (14), and similar to 
the discussion in Section |2.1[ we have 


z (t) =Az ( t ) + B £ ( t ), 
£ (o) (t) =Cz{t) + Di{t). 


A, B , C and D are real matrices 


where 


A = R-'y(X-h) 

B = [ ~^XS 2l -y/Rl % } 

C = ^7 S 12 X 

D = (Sii + S12XS21) [ h C> 4 x 8 ] 


( 21 ) 


( 22 ) 


R = 


X 

Su 

5 21 
S 12 

5 22 


7+ K 

2 


0 

7+k 

2 

0 

_ 6 

2 

0 

0 

0 

0 


(Is — S 22 ) 1 , 

[ I 4 ©4x8 ] s 

[ 0$x4 h } S 

[ h ©4x8 ] s 

[ 08x4 h ] S 


2 

0 

2 

0 

0 

0 

0 

0 


h 

Osx4 

h 

08x4 

O 4 X 8 

h 

O 4 X 8 

h 


0 

_e 

2 

0 

7+k 

2 

0 

0 

0 

0 


0 

0 

0 

0 

7+« 

2 

0 

2 

0 


0 

0 

0 

0 

0 

7+K 

2 

0 

_ € 

2 


0 

0 

0 

0 

€ 

2 

0 

7 +K 
2 

0 


0 

0 

0 

0 

0 

_ € 

2 

0 

7+K 


(23) 


More details of how to obtain (21) are given in Appendix 1. Based on (§, @ and 
we have 


H 

H 2 


D - CA~ l B 
[1 0 1 0 ]H, 
[0 1 0 -1 ] H, 


(24) 





















and the two-mode squeezing spectra are 


where 


V (0) = U+ (0) + V.- (0) = Tr [H*H\ + H* 2 H 2 ] = TV [ H*M h2 H} , 


Mi t2 


10 10 

0 10-1 
10 10 

0-101 


(25) 


(26) 


Since V (0) is parametrized by the matrix S or S, we define V(0; S) as the value of V (0) 
for a fixed value of S, and V ( 0 ; S) as the value of V ( 0 ) for a fixed value of S. 


4 Optimization of S 


Using the same parameters (total pump power and damping rates) as that of the coherent- 
feedback dual NOPA system described in [7], for each NOPA we set the constant relating 
to the amplitude of the pump beam e = 0 . 47 r and the damping rate 7 = 7 r , where 
7 r = 7.2 x 10' Hz is the reference value for the transmissivity rate of the mirrors. We 
consider the system in the ideal case, where there are no losses (k = 0). Following Section 
[3] we compute the two-mode squeezing of the dual-NOPA coherent feedback network to 
be 10 log 10 V (0; S cfb ) = -26.235 dB. 


In this section, we aim to find a complex unitary matrix S) m at which the cost function 
V (0; S) is locally minimized. We will numerically solve this optimization problem using 
the method of modified steepest descent on a Stiefel manifold, which reformulates the 
problem with a unitary constraint as an unconstrained problem on a Stiefel manifold. The 


Stiefel manifold in our problem is the set St( 6,6) = jS 1 e C 6x6 : S*S — /j. The modified 
steepest descent on a Stiefel manifold method employs the first-order derivative of the cost 
function, more details about this algorithm can be found in [[15] . 

For any square matrix Y such that / — Y is invertible, we have (/ — Y)~ l = (/ — 
U) _1 (J + Y — Y) = / + (/ — Y) _1 Y. Based on the above fact and equations (40), (22) and 
(24), we expand H(S + AS) as H(S) + H(AS ) + 0(AS' 2 ), where 0(AS' 2 ) denotes terms 
that depend on terms that are products containing at least two AS. 

By using (25), we have 


V^S + AS) = Tr[H(S + AS)*M lj2 H(S + AS)} 

= V (0; S ) + Tr[H(AS)*M h2 H(S) + H(S)*M 1>2 H(AS)\ + 0(|| AS|| 2 ),(27) 


where 0(||AS'|| 2 ) denotes that the function 0(||AS'|| 2 ) satisfies 1 < c for some 

positive constant c for all ||AS'|| > 0 sufficiently small. Since H(AS ) and H(S) are real 
matrices at cu = 0 , Mi _ 2 is a real symmetric matrix, and a matrix and its transpose have 
the same trace, we get 


U(0;5 +A5) = U(0;S') + 2Ti[MH(AS)] + 0(||A5|| 2 ), (28) 
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where M = H(S)*Mi j2 ■ Based on (19) and the property that trace is invariant under cyclic 
permutations, we have 


V(0;S + AS) = R(0; S) + ReTV[AS*Do] + 0(||AS'|| 2 ), 


(29) 


where 


N = 


2 K*N*K 

( 

I 4 

V 

o 8x4 


-V7 


(30) 


[ h 04X8 ] M + °^ S XS 2 1 [ h 0 4 xs ] M 

0 4x 8 

h 


xa~ 1 bm\ ([ h o 4x8 ] +S 12 X [ 0 8x4 Is ] 


+^CA~ 1 X [ 0 8x4 Is]). 


(31) 


When AS approaches the zero matrix 06,6, D§ is the directional derivative of 11(0; S ) at 
S in the direction AS H5l. 


By applying the modified steepest descent on Stiefcl manifold method, we use the 
following steps to find a matrix S) m at which the dual-NOPA system is stable and V (0; S) 
is locally minimized. The system is stable if and only if the matrix A in equation (21) 


is Hurwitz, that is, real parts of all the eigenvalues of A are negative. Denote the vector 
of all the mode operators of our system by z(t) = [a 4 (f), b\(t), a 2 (t), b 2 {t)] T and define the 
intra-cavity photon number operator as n(t) = z(t)*z(t). If the system is stable, we have 
lim^oo(£(£)) = 0 [9]. Moreover, based on the quantum Ito rules [T6], Section 2.5 in mi 
verifies that lim^ 00 (n(t)) is bounded. Therefore, a stable system physically means that 
the mean total number of intra-cavity photons in the system would not keep increasing as 
time t approaches infinity. Note that, based on (22) and (24), to get V(0;S'), (Is — S 22 ) 
and A must be invertible. The algorithm is initiated at the point S c fb given by (17), the 
static passive linear matrix denoting the dual-NOPA coherent feedback network shown in 

Fig. [T} 

Step 1. Start from S = S c fb, choose step size p — 1. 

Step 2. 

and compute the descent direction Z = SD*-S — D 

Step 3. Calculate (Z,Z) = Tr[Z*(I — \SS*)Z\. If a J(Z, Z) is smaller than 10 
and set Si m = S. 

Step 4■ Calculate Si = tt(S + 2pZ) (Note: if the singular value decomposition (SVD) of 
a matrix X is X = I/SB*, then tt(X) = UIV*). Calculate A and (J 8 — S 22 ) corresponding 
to S\ based on ( J 22 ] ) and ([24). If max(Re(eig(A))) > 0 or det(A) = 0 or det(/ 8 — S 22 ) = 0 
go to step 5. Otherwise, if V(S) — B(S'i) > p(Z, Z), set 2 p —> p, and go back to Step 4. 

Step 5. Calculate S 2 = tt(S + pZ). Compute A and (I 8 — S 22 ) corresponding to S 2 . 
While max(Re(eig(A))) > 0 or det(A) = 0 or det(/ 8 — ^ 22 ) = 0, set \p -A p, go back to 
step 5. If V(S) — V(S 2 ) < \p{Z, Z) , set \p -A p, and go back to Step 5. 

Step 6. Set S = S 2 and repeat Step 2. 


Calculate Dg , the directional derivative of R(0;S), by using equation (30), 


4-3 


stop 
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Thus, we get S[ m and the corresponding gradient G = D§ shown as (44) and (|45j) in 
Appendix 2, the operator norm of G is || G ||= 6.497 x 10~ 4 , the norm of a tangent direction 
Z is \J(Z, Z) = 9.112 x 10” 4 and EPR entanglement is locally maximized for a locally 
minimum value of V (0; S). The locally minimized two-mode squeezing is 10 log 10 V (S) m ) = 
10 log 10 (4.1824x 10 -8 ) = —73.786 dB, which is about 47.551 dB less than the value reported 
in [7] for the dual NOPA coherent feedback system in Fig. [I] for the same values of 
the parameters of the NOPAs. Also, it is checked by equations ((8|) , Q and ( f24| ) that, 


5 Decomposition of unitary matrices 


As introduced in Section [TJ the 6x6 matrix S) m denotes a network constructed from static 
passive linear optical devices. To find the specific physical configuration of the network, 
we first employ the approach in [T8j to decompose S) m as the product of 15 two-level 6x6 
unitary matrices, that is, Si m = II^SA Here, a two-level n x n unitary matrix refers 
to a special type of unitary matrix that has a unitary 2x2 principal submatrix and the 
remaining matrix elements are the same as those of the n x n identity matrix. Each of 
these two-level unitary matrices has determinant with modulus 1. The reason we use this 
method here is that any two-level unitary matrix is isomorphic to the set of 2 x 2 unitary 
matrices, which represent the transformation performed in the Heisenberg picture by static 
linear optical devices, such as beam splitters and phase shifters. 


The decomposition of Si m does not give a unique group of two-level unitary matri¬ 
ces. Types of two-level unitary matrices in a group are determined by a vector P = 
(p i, p 2 , • • • , Pe) j where entries correspond to a permutation of (1,2, ••• ,6). A two-level 
matrix is named P-unitary matrix of type k if row and column indexes of its principal 
submatrix are pk and Pk+i■ By setting P = (6, 5,4, 3, 2,1) and using the Matlab program 
pub.m developed by [18], we find a group of two-level unitary matrices Sk (k G [1,13]) 
shown as (46) in Appendix 2. 

A unitary matrix representing a beamsplitter has the form [19] 


a P 

—/3 a 


(32) 


where transmission rate a and reflection rate f3 are real numbers satisfying |a| J + \/3 \ 2 = 
1. Thus, Sk (k € {3,5,6,8,12,13}) represents transformation by a beamsplitter, with 
parameters 


a 3 = 

-ai 3 = 0.9999632197, 


Ot 8 = 

-a 5 = 0.0084711563, 


Oil2 = 

-a 6 = 0.0123787627, 


Pk = 

Jl-al ke{ 3,5,6,8,12,13}. 

(33) 


The configuration of the whole network is shown in Fig. |4} The network requires high 
accuracy of the value of the parameter ct*,. To achieve 10 log 10 V(Si m ) = —73.786 dB, we 
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need to keep at least six decimal places for ctk- However, with a lower accuracy of less than 
six digits but more than one digit, we still get better two-mode squeezing than that of 
the dual-NOPA coherent feedback network. For example, by rounding off to two decimal 
places, that is, ct 3 = — aq 3 — 1 , a 8 — a 12 = —a 5 = — a 6 = 0 . 01 , with /3 fc = ^1 — a\ 
(k e {3, 5, 6 , 8,12,13}) as before, the two mode squeezing of the optimized network is 
—36.546 dB. When accuracy is less than two decimal digits, the network becomes exactly 
the dual-NOPA coherent feedback network. 



Figure 4: Physical configuration of the optimized dual-NOPA network. The system con¬ 
tains NOPAl, NOPA2 and a static passive linear network which consists of six beam 
splitters Sk, k G (3, 5, 6 , 8,12,13} (denoted by the red lines) with the parameters as 
shown in equation (33). For each NOPA, n = 0, e = 0 . 47 r , 7 = 7 r , where 7 r = 7.2 x 10'. 
The black lines are mirrors that are fully reflecting. Inputs fields 1 and £j n) 2 are in the 
vacuum state and EPR entanglement is generated between the two outputs £out,i and ^ 0 ut .,2 
as discussed in Section U 


6 Conclusion 

This paper has studied the optimization of EPR entanglement in terms of maximising the 
two-mode squeezing generated by a quantum network that contains two NOPAs connected 
by a static passive linear network. The transformation implemented by the passive network 
is a 6 x 6 complex unitary matrix. By employing the modified steepest descent on Stiefel 
manifold method, we have found the passive network 5) m at which the two-mode squeez¬ 
ing function V (0; S) used to evaluate the EPR entanglement is approximately minimized 
locally. It is shown that E + (0;S) m ) ~ H_(0;S) m ) and H(0;S) m ) = 4.1824 x 1CT 8 , which 
approximates the ideal case of infinite squeezing where V+(0; S) = R_(0;S') = 0. Also, 
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with the same values of the parameters e and 7 and without considering losses (k — 0 ), 
the optimized network improves the EPR entanglement by a significant reduction of 47.551 
dB in the two-mode squeezing compared to the one of the dual-NOPA coherent feedback 
network of Fig. [l] studied in [7|. Finally, by decomposing S) m into a product of two-level 
unitary matrices, we have found the physical set up of the optimized network as shown in 
Fig. i The network requires that beamsplitters have highly accurate realization of the 
transmission rates a*,, k E {3, 5, 6,8,12,13} with an accuracy of at least six decimal places 
to achieve the two-mode squeezing 10 log 10 V (0; S\ m ) = —73.786 dB. If accuracy is less than 
six and more than one decimal digits, the two-mode squeezing is lower than —73.786 dB 
but better than that of the dual-NOPA coherent feedback network in Fig. [lj If accuracy 
of is zero or one decimal digit, the network coincides with the dual-NOPA coherent 
feedback network in Fig. [l] 


APPENDIX 1 

For convenience, define the following quadratures 

* =K, < b\, 6?, al,a p 2 , b\, 

p. = \p1 fP pi pP pi pP £9 tP 1 T 

s«n [>>in,a,l ’ >m,a, 1 ’ 1 ’ 5 > in,a,2 ’ >m,a,2 ’ 5 2-1 ’ 

p += \& pP £P p* pP pQ pP ] T 

£,out [>>out,a,li ^>out,a, 1’ ^>out,b, 1’ ( *out,b,l'> ^>out,a, 2’ >cmf,a,2’ ^>out,b, 2’ ^>out,b,2i ’ 

P, = \p1 fP fl fP pi PP Pi PP ] T 

Stoss i^>loss,a,l 5 Moss,a,1 ’ Moss,6,1 ’ Moss,6,1 ’ ^ loss,a,2 ’ Moss,a,2 ’ ^ loss,b,2 ’ Moss,6,2-1 ’ 

AO = rt9 pP pi pP ]T do) = \ti pp Pi cP 1T 
^ l.Mn,l ’ >m,l ’ Mn,2’Mn,2J ’ ^ LSont,lJ ^oit£,l’ ^out,2'> ^out,2i > 

Cl=K ( °\Cm] T , 6 = [£ W ,£out] T , 

C 7 = [£m,6oss] T , £ = 6 oss] T - 


According to (13) and (14), we have 


i (f) = (t) + [ -y^y/s —\/Ki8 ] £' (0, 

iout (t) = yf^Z (t) + iin (t) , 


where R is as shown in (23) 
Define 


6 CO = ^ 6 ( 0 , 


(34) 


(35) 

(36) 


(37) 


in which the real unitary matrix S is the quadrature form of S, satisfying the following 
relations 


S T S = SS T = I 12 , 
9/ q t — t 

u u n u — *- , ni 


where 


Jn — h 


0 1 

-1 0 


(38) 


(39) 
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Let 


S = 


Sn s 12 

*521 *522 


(40) 


with 


We have 


Su = 

h 04x8 

5 

1 1 1 

X 

1—1 00 

0 

5 

S21 = 

X 

00 

0 


0 ^ 

£ ^ 
4 ^ 

5 

S12 = 

O 

4 ^ 

X 

00 


1 1 

00 

O 

5 

S22 = 

X 

00 

0 

5 

1 1 

00 

0 

(41) 


£ (o) = Snt® + S 12 Ut, 
£in = *5 2 1*£ W + S22^out- 


Based on (35), (36), (42) and (43), we obtain (21). 


(42) 

(43) 


APPENDIX 2 


Si m and its corresponding gradient G = Dg are as shown in (44) and (45), respectively. 
The group of two-level unitary matrices Sk (k G [1,13]) as the decomposition of S) m in 


'5/m — 


- 0.012305658659326 

- 0.000000000000071 

0.999887502042829 

0.000000000000051 

0.008576364236157 

-0.000000000000001 


0.000000000000071 

- 0.012305658659326 

0.000000000000001 

0.008576364236158 

- 0.000000000000050 

0.999887502042829 


0.008576364236157 

- 0.000000000000109 

- 0.008471156255372 

0.000000000000085 

0.999927340104363 

0.000000000000047 


- 0.000000000000142 

0.999887502042830 

0.000000000000069 

- 0.008471156255372 

0.000000000000111 

0.012378318554963 


0.999887502042829 

0.000000000000142 

0.012378318554964 

-0.000000000000112 

- 0.008471156255372 

- 0.000000000000069 


0.000000000000110 

0.008576364236157 

- 0.000000000000048 

0.999927340104363 

- 0.000000000000085 

- 0.008471156255373 


( 44 ) 
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G = 1(T 3 


0.409017637139186 

0.000013384124031 

0.007586926959628 

-0.000009579672002 

-0.292783525876746 

0.000000246886333 


-0.000013389581000 

0.409017646234133 

-0.000000246987169 

-0.292783532386752 

0.000009583578405 

0.007586927128612 


-0.292752071754183 

-0.092604224236176 

-0.005430300582194 

0.066288073024820 

0.209558161024828 

-0.001717728084123 


-0.131726307590905 

0.007631674074983 

-0.002443410193280 

-0.005462914965842 

0.094292493228142 

0.000141561463910 


0.007542180013146 

0.131726307489625 

0.000139901426214 

-0.094292493155211 

-0.005398853213792 

0.002443410191339 


0.092604228142405 

-0.292814986508976 

0.001717728155075 

0.209603196744500 

-0.066288075820747 

-0.005431467597355 


(45) 


Section [5] is 

S\ = S15 = Iq, 

'10 0 0 

0 0 10 


0 0 0 0 

0 0 0 

' 1 0 0 0 0 0 ' 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 0 1 0 

0 0 0 -1 0 0 

0 0 0 0 0 1 _ 

chq / 3 g 0 0 0 0 

—/?6 Ctg 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0’ 

0 0 0 0 1 0 

0 0 0 0 0 1 _ 

-1 0 0000' 

0 -1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0 ’ 

0 0 0 0 1 0 

0 0 0 0 0 1 _ 

«13 /?13 0 0 0 0 

—/3i3 «i3 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0' 

0 0 0 0 1 0 

0 0 0 0 0 1 


0 

0 ' 


1 

0 

0 

0 

0 

0 


0 

0 


0 

1 

0 

0 

0 

0 


0 

0 

0 

0 

, £3 = 

0 

0 

0 

0 

«3 

-/3.3 

/3s 0 

a 3 0 

0 

0 


1 

0 


0 

0 

0 

0 

1 

0 


0 

1 _ 


0 

0 

0 

0 

0 

1 





' 1 

0 

0 0 

0 

0 

- 




0 

1 

0 0 

0 

0 





0 

0 

1 0 

0 

0 




Sfj = 

0 

0 

0 1 

0 

0 





0 

0 

0 0 

a 5 

& 





_ 0 

0 

0 0 

5 

a 5 

. 




' 1 

0 

0 

0 

0 

0 

" 




0 

1 

0 

0 

0 

0 




S 8 = 

0 

0 

0 

0 

«8 
—Ps 

P "Cfc 

00 00 

O O 

0 

0 





0 

0 

0 

0 

1 

0 





_ 0 

0 

0 

0 

0 

1 





' 1 

0 

0 

0 

0 

0 




0 

1 

0 

0 

0 

0 



Sl 2 = 

0 

0 

0 

0 

a 12 
~Pl 2 

P 12 
«12 

0 

0 

0 

0 




0 

0 

0 

0 

1 

0 




_ 0 

0 

0 

0 

0 

1 
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